#include <iostream>
#include <vector>
#include <cmath>
#include <fstream>
#include <string>
#include <iomanip>  
#include <cblas.h>
#include <lapacke.h>
#include "Spline.hpp"  

// Define the function f(x) = 1 / (1 + 25 * x^2)
double f(double x) {
    return 1.0 / (1.0 + 25 * x * x);
}

// Calculate the error between the exact function and the interpolated spline
double computeError(const PPForm& spline, const std::vector<double>& x_vals, const std::vector<double>& y_vals) {
    int n = x_vals.size();
    double maxError = 0.0;

    // Compute error at midpoints of each interval
    for (int i = 0; i < n-1; ++i) {
        double midpoint = (x_vals[i] + x_vals[i+1]) / 2.0;
        double exact_val = f(midpoint);
        double interp_val = spline.getSplineValue(midpoint);
        maxError = std::max(maxError, std::abs(exact_val - interp_val));
    }

    return maxError;
}

// Run the test for different values of N and store results
void runTest(const std::vector<int>& N_values, const std::string& boundary_type, double first_der_a = 25/26.0/26.0, double first_der_b = -25/26.0/26.0) {
    std::cout << "N, Max Error\n";

    for (int N : N_values) {
        // Generate evenly spaced nodes and calculate function values
        std::vector<double> x_vals(N+1), y_vals(N+1);
        
        // Write data to a file for each N
        std::ofstream file("spline_A_" + std::to_string(N) + ".csv");
        for (int i = 0; i <= N; ++i) {
            x_vals[i] = -1.0 + 2.0 * i / (N);  // Evenly spaced nodes from -1 to 1
            y_vals[i] = f(x_vals[i]);
        }
        // Create the PPForm spline interpolation object
        PPForm spline(x_vals, y_vals, boundary_type, first_der_a, first_der_b);
        for (double x = -1; x <= 1; x += 0.01) {
            double y = spline.getSplineValue(x);  // 通过样条插值计算 y 值
            file << x << "," << y << "\n";  // 将结果写入文件
        }
        // Compute the maximum error at midpoints
        double maxError = computeError(spline, x_vals, y_vals);
        file.close();
        std::cout << "N = " << N << ", Max Error = " << maxError << "\n";
    }
    std::cout << "Test results saved to spline_errors.csv" << std::endl;
}

int main() {
    // Define the different N values for testing
    std::vector<int> N_values = {6, 11, 21, 41, 81};

    // Boundary conditions types: "natural", "periodic", or "complete"
    std::string boundary_type = "complete";  // You can change this to "periodic" or "complete" based on your needs

    // Run the test for different boundary conditions
    runTest(N_values, boundary_type);

    return 0;
}
